/*
 * real_li.c
 *
 * Computation of the logarithmic integral for real values.
 *      Author: Kevin Marshall Stueve; August 2014
 * Copyright (C) 2014 Kevin Marshall Stueve
 *
 * This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 *
 */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "real_li.h"

/*
 * A test of the logarithmic integral code.  The agreement between two implementations verifies the code.
 *     Authors: Fredrik Johansson; Summer 2009; python version
 *              Kevin Marshall Stueve; August 2014; c version
 */
void real_li_test() {
	printf("Real logarithmic integral test on y=10^n for n between 1 and 11, inclusive:\n");
	for (int n=1;n<12;n+=1) {
	    long double y = pow(10,n);
	   	printf("n:%u,y:%Lf\n",n,y);
	   	printf("li_approx(y)-li2: %Lf, li(y)-li2: %Lf\n", li_approx(y)-li2, li(y)-li2);
	   	printf("\n");
	}

}

/**
 * Compute the absolute value of x.
 *     Author: Kevin Marshall Stueve; August 2014
 */
double double_abs(double x){
	if (x>=0) {
		return x;
	}
	return -x;
}
/*
 * The logarithmic integral function, li(x), is not to be confused with the offset logarithmic integral function, Li(x).
 * They are defined as follows:
 * li(x)=\int ^x_0 \frac{dt}{log(t))}
 * Li(x)=\int ^x_2 \frac{dt}{log(t))}
 * Li(x) = li(x)-li(2)
 * An efficient way to calculate li(x) is with the exponential integral function:
 * li(x)=Ei(log(x))
 * The series rerpesention for Ei can be found on Wikipedia:
 *  li(e^u) = Ei(u) = \gamma + \log |u| + \sum_{n=1}^\infty {u^{n}\over n \cdot n!}
 *
 * Code for computation of the logarithmic integral based on the exponential integral was provided by Fredrik Johansson; Summer 2009.
 * The code is below.
*/
#define gamma 0.57721566490153287
#define ei_eps 1e-17
#define asymp_cutoff 50//50
//#define li2 1.04516378011749278 //li(2) (defined in shared.h)
/*
 * Computes the exponential integral which is used to compute the logarithmic integral.
 *     Authors: Fredrik Johansson; Summer 2009; python version
 *              Kevin Marshall Stueve; August 2014; c version
 */
long double ei(long double x) {
    long double xexp=exp(x);
    long double s;
	long double t;
	int k;
	if (abs(x) < asymp_cutoff){
		s = gamma + log(double_abs(x));
		t = 1.0;
		k = 1;
		while (double_abs(t) > ei_eps) {
			t *= x;
			t /= k;
			s += t / k;
			k += 1;
		}
	} else {
	    s = 1;
	    t = 1.0;
		k = 1;
		long double r = 1. / x;
	    while (double_abs(t) > ei_eps){
	    	t *= k;
	    	t *= r;
	    	s += t;
	    	k += 1;
	    	s = s * r * xexp;
	    }
	}
    return s;
}
/*
 * The logarithmic integral function.
 *     Authors: Fredrik Johansson; Summer 2009; python version
 *              Kevin Marshall Stueve; August 2014; c version
 */
long double li(long double x){
	return ei(log(x));
}
/*
 * An asymptotic series for the logarithmic integral that is a close approximation
 *     Authors: Fredrik Johansson; Summer 2009; python version
 *              Kevin Marshall Stueve; August 2014; c version
 */
long double li_approx(long double y){
	long double x = log(y);
	long double t=y;
    long double s=t;
	long double r = 1. / x;
	int k = 1;
    long double tprev = t;
    while (t <= tprev) {
        tprev = t;
        t *= k;
        t *= r;
        s += t;
        k += 1;
    }
    return s*r;
}
/*
 * The "Offset" logarithmic integral, where the integration starts at 2 instead of 0.
 *     Authors: Fredrik Johansson; Summer 2009 python version
 *              Kevin Marshall Stueve; August 2014; c version
 */
long double Li2(long double x){
	return li(x) - li2;
}
